/*
**--------------------------------------------------------------
**Estimation of a spatio-temporal autoregressive (STAR) model
**
** Reference: 
**     Dube and Legros (2014) - 
**				Letters in Spatial and Resource Sciences
**	   Thanos et al. (2016) - 
**				Regional Science & Urban Economics
**
** Written by
**  Jean Dube (U. Laval) 
**  Diego Legros (U. Bourgogne)
**	Nicolas Devaux (UQAR)
**  Sotirios Thanos (U. Manchester)
**--------------------------------------------------------------
*/

**--------------------------------------------------
**Preparing Stata for estimation
**--------------------------------------------------
clear mata
clear all

**Loading the data base
use "Transaction-STExample.dta"

**Generating temporal (continuous) variable
quietly generate date = 365*(year - 1990) + 31*(month - 1) + day
sort date   /*Chronologically order data*/

**Generating coordinates in km (instead of meters)
** (eliminating scale problem for building weights)
quietly generate xc = XMTM/1000
quietly generate yc = YMTM/1000

**Generating dependent variable
quietly generate log_price    = log(saleprice)

/*Generating independent variables*/

**Generating group of independent dummy variables
tab quality, gen(QIndex)
tab location, gen(Sfixed)

**Generating time fixed effects (annual dummy variables) 
quietly generate d91 = (year==1991)
quietly generate d92 = (year==1992)
quietly generate d93 = (year==1993)
quietly generate d94 = (year==1994)
quietly generate d95 = (year==1995)

**Applying mathematical transformation on independent variables
quietly generate log_livearea = log(livearea)
quietly generate log_age = log(age)

**--------------------------------------------------
**Build spatio-temporally lagged variables
**--------------------------------------------------
quietly generate WTprice = .

mata 
	/*Import the main vectors used for calculations*/
	N  = st_nobs()
	yn = st_data(.,"log_price")    /*Vector of dependent variables*/
	XC = st_data(.,"xc")           /*Vector of X coordinates*/
	YC = st_data(.,"yc")           /*Vector of Y coordinates*/
	TC = st_data(.,"date")         /*Vector of temporal coordinates*/
	/*Generate the spatial and temporal matrices*/
	XX = (XC*J(N,1,1)') - (J(N,1,1)*XC')  /*Distance between X coordinates*/
	YY = (YC*J(N,1,1)') - (J(N,1,1)*YC')  /*Distance between Y coordinates*/
	ZZ = (TC*J(N,1,1)') - (J(N,1,1)*TC')  /*Eucledian distance between transactions*/
	DD = sqrt((XX:*XX) + (YY:*YY))        /*Distance matrix*/
	SO = exp(-DD)-I(N)         /*Spatial weight matrix (negative exponential transformation)*/
	TT = exp(-abs(ZZ))         /*Temporal weight matrix (negative exponential transformation)*/
	/*Impose restrictions on spatial and temporal relations*/
	DK = ((DD*J(N,1,1)):/N)*J(N,1,1)'      /*Estimate & attribute to each observtion the mean distance of all observations*/
	T1 = (ZZ:<=65):*(ZZ:>25)   /*Define the temporal "past" window (in this case between 25 and 65 days)*/
	/*Generating spatial and spatio-temporal weights*/ 
	ST = (DD:<=DK):*SO         /*Truncated spatial weights matrix*/
	TU = T1:*TT                /*Temporal weights within temporal "past" window*/
	Wb = ST:*TU                /*Lower triangular spatio-temporal weight matrix*/
	/*Standardizing spatio-temporal weight matrix*/
	SL  = (Wb*J(N,1,1))*J(N,1,1)' /*Row sums*/
	WSb = Wb:/(SL + (SL:==0))     /*Row standardization*/
	/*Generating the spatiotemporally dynamic depedent variable*/
	Wyn = WSb*yn			     /*Creating spatio-temporally lagged dependent variable (exogenous variable)*/
	st_store(.,("WTprice"),(Wyn))
end

**Eliminating observations having no spatio-temporal neighbours
drop if WTprice==0
clear mata

**--------------------------------------------------
**Building spatial block-diagonal weights matrix
**--------------------------------------------------
mata 
	/*Import the main vectors used for calculations*/
	N  = st_nobs()
	yn = st_data(.,"log_price")    /*Vector of dependent variables*/
	XC = st_data(.,"xc")           /*Vector of X coordinates*/
	YC = st_data(.,"yc")           /*Vector of Y coordinates*/
	TC = st_data(.,"date")         /*Vector of temporal coordinates*/
	/*Generate the spatial and temporal matrices*/
	XX = (XC*J(N,1,1)') - (J(N,1,1)*XC')  /*Distance between X coordinates*/
	YY = (YC*J(N,1,1)') - (J(N,1,1)*YC')  /*Distance abetween Y coordinates*/
	ZZ = (TC*J(N,1,1)') - (J(N,1,1)*TC')  /*Eucledian distance between transactions*/
	DD = sqrt((XX:*XX) + (YY:*YY))        /*Distance matrix*/
	SO = exp(-DD)-I(N)         /*Spatial weight matrix (negative exponential transformation)*/
	TT = exp(-abs(ZZ))         /*Temporal weight matrix (negative exponential transformation)*/
	/*Impose restrictions on spatial and temporal relations*/
	DK = ((DD*J(N,1,1)):/N)*J(N,1,1)'      /*Estimate & attribute to each observtion the mean distance of all observations*/
	TO = (ZZ:<=25):*(ZZ:>=-15) /*Define the "present" temporal window or contemporaneous period*/
	/*Creating contemporaneous spatial weight matrix*/
	TO = TO:*TT
	Sb = SO:*TO
	/*Row-standardizing the contemporaneous weight matrix*/
	SL  = (Sb*J(N,1,1))*J(N,1,1)' /*Row sums*/
	SSb = Sb:/SL                  /*Row standardization*/
	/*Exporting the contemporaneous weight matrix*/
	st_matrix("Sb",Sb)
end

**--------------------------------------------------
**Importing spatial weights matrix in Stata
**--------------------------------------------------
spmat putmatrix WST Sb, normalize(row) id(ID) replace
spmat note WST : "Spatial - block diagonal - row-standardized weights matrix"
spmat graph WST, blocks(10)  /*Graphing the spatial relations of the weights matrix*/

**--------------------------------------------------
**Generating macro for independent variables
**--------------------------------------------------

global varX log_livearea log_age basement garage panoramicview QIndex2 QIndex3 bathroom
global varT d91-d95
global varS Sfixed1 Sfixed2 Sfixed4 Sfixed5

**--------------------------------------------------
**Estimating OLS model
**--------------------------------------------------

reg log_price $varX $varT 
predict error, resid

/*Graphing residuals - vizual detection of spatial autocorrelation*/
#delimit ;
twoway 
	(scatter YMTM XMTM if error>0, mcolor(red) msize(medium) msymbol(lgx)) 
	(scatter YMTM XMTM if error<0, mcolor(blue) msize(medsmall) msymbol(plus)), 
	title(Location of the prediction error) 
	subtitle(Red: positive error - Blue: negative error) 
	legend(off);
#delimit cr

**--------------------------------------------------
**Estimating STAR models
**--------------------------------------------------

/*STAR model - spatial unidirectional effect*/
reg log_price $varX $varT WTprice 

/*SAR model - spatial multidirectional effect*/
spreg gs2sls log_price $varX $varT, id(ID) dlmat(WST)

/*STAR model - spatial unidirectional & multidirectional effects*/
spreg gs2sls log_price $varX $varT WTprice, id(ID) dlmat(WST)
